HW_1
Dmitrii Belousov
2024-11-16
Data Import
data <- readRDS("very_low_birthweight.RDS")Task 1
Загрузите датасет very_low_birthweight.RDS (лежит в папке домашнего
задания). Это данные о 671 младенце с очень низкой массой тела (<1600
грамм), собранные в Duke University Medical Center доктором Майклом О’Ши
c 1981 по 1987 г.
Описание переменных см. здесь. Переменными исхода являются колонки
‘dead’, а также время от рождения до смерти или выписки (выводятся из
‘birth’ и ‘exit’ 7 пациентов были выписаны до рождения). Сделайте копию
датасета, в которой удалите колонки с количеством пропусков больше 100,
а затем удалите все строки с пропусками.
data %>%
head()## birth exit hospstay lowph pltct race bwt gest inout twn lol
## 1 81.511 81.604 34 NA 100 white 1250 35 born at Duke 0 NA
## 2 81.514 81.539 9 7.250000 244 white 1370 32 born at Duke 0 NA
## 3 81.552 81.552 -2 7.059998 114 black 620 23 born at Duke 0 NA
## 4 81.558 81.667 40 7.250000 182 black 1480 32 born at Duke 0 NA
## 5 81.593 81.599 2 6.969997 54 black 925 28 born at Duke 0 NA
## 6 81.602 81.771 62 7.189999 NA white 940 28 born at Duke 0 NA
## magsulf meth toc delivery apg1 vent pneumo pda cld pvh ivh ipe
## 1 NA 0 0 abdominal 8 0 0 0 0 <NA> <NA> <NA>
## 2 NA 1 0 abdominal 7 0 0 0 0 <NA> <NA> <NA>
## 3 NA 0 1 vaginal 1 1 0 0 NA <NA> <NA> <NA>
## 4 NA 1 0 vaginal 8 0 0 0 0 <NA> <NA> <NA>
## 5 NA 0 0 abdominal 5 1 1 0 0 definite definite <NA>
## 6 NA 1 0 abdominal 8 1 0 0 0 absent absent absent
## year sex dead
## 1 81.51196 female 0
## 2 81.51471 female 0
## 3 81.55304 female 1
## 4 81.55847 male 0
## 5 81.59406 female 1
## 6 81.60229 female 0cols_to_select <- data.frame(NA_count = colSums(is.na(data))) %>% filter(NA_count <= 100) %>% row.names()
data %>%
select(all_of(cols_to_select)) %>%
na.omit() %>%
write_csv("filtered_dataset.csv")Task 2
Постройте графики плотности распределения для числовых переменных. Удалите выбросы, если таковые имеются. Преобразуйте категориальные переменные в факторы. Для любых двух числовых переменных раскрасьте график по переменной ‘inout’.
data <- data %>%
mutate(across(where(is.character), as.numeric),
across(where(is.factor), as.factor),
across(where(~ is.numeric(.) & n_distinct(., na.rm = TRUE) == 2), as.factor))
data %>%
select_if(is.numeric) %>%
pivot_longer(cols = everything(), names_to = "variable", values_to = "value") %>%
ggplot(aes(x = value)) +
geom_histogram(bins = 10, fill = "skyblue", color = "black") +
facet_wrap(~variable, scales = "free_x") +
theme_minimal() +
theme(axis.text.x = element_text(size=7))## Warning: Removed 657 rows containing non-finite outside the scale range
## (`stat_bin()`).
dfSummary(data)## Data Frame Summary
## data
## Dimensions: 671 x 26
## Duplicates: 0
##
## --------------------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ----------- ---------------------------- --------------------- --------------------- ---------- ---------
## 1 birth Mean (sd) : 84.8 (1.6) 542 distinct values : : 650 21
## [numeric] min < med < max: : : : : . : (96.9%) (3.1%)
## 81.5 < 84.9 < 87.5 : : : : : : : : :
## IQR (CV) : 2.6 (0) : : : : : : : : : :
## : : : : : : : : : :
##
## 2 exit Mean (sd) : 84.8 (1.8) 539 distinct values : 640 31
## [numeric] min < med < max: : . (95.4%) (4.6%)
## 68.5 < 85 < 96.9 : :
## IQR (CV) : 2.6 (0) : :
## : : :
##
## 3 hospstay Mean (sd) : 40.4 (304.8) 153 distinct values : 640 31
## [integer] min < med < max: : (95.4%) (4.6%)
## -6574 < 37 < 3668 :
## IQR (CV) : 46 (7.6) :
## :
##
## 4 lowph Mean (sd) : 7.2 (0.1) 74 distinct values : . 609 62
## [numeric] min < med < max: : : (90.8%) (9.2%)
## 6.5 < 7.2 < 7.5 . : :
## IQR (CV) : 0.2 (0) . : : : :
## . : : : : : .
##
## 5 pltct Mean (sd) : 201.6 (80.5) 266 distinct values : 601 70
## [integer] min < med < max: . : . (89.6%) (10.4%)
## 16 < 202 < 571 . : : :
## IQR (CV) : 109 (0.4) : : : : .
## : : : : : : .
##
## 6 race 1. black 369 (57.1%) IIIIIIIIIII 646 25
## [factor] 2. native American 16 ( 2.5%) (96.3%) (3.7%)
## 3. oriental 4 ( 0.6%)
## 4. white 257 (39.8%) IIIIIII
##
## 7 bwt Mean (sd) : 1093.9 (265.2) 142 distinct values . : 669 2
## [integer] min < med < max: . : : : : (99.7%) (0.3%)
## 400 < 1120 < 1580 : : : : : :
## IQR (CV) : 410 (0.2) . : : : : : : : :
## . : : : : : : : : :
##
## 8 gest Mean (sd) : 28.9 (2.5) 24 distinct values : : 667 4
## [numeric] min < med < max: : : : (99.4%) (0.6%)
## 22 < 29 < 40 . : : :
## IQR (CV) : 4 (0.1) : : : :
## . : : : : :
##
## 9 inout 1. born at Duke 547 (81.9%) IIIIIIIIIIIIIIII 668 3
## [factor] 2. transported 121 (18.1%) III (99.6%) (0.4%)
##
## 10 twn 1. 0 516 (79.3%) IIIIIIIIIIIIIII 651 20
## [factor] 2. 1 135 (20.7%) IIII (97.0%) (3.0%)
##
## 11 lol Mean (sd) : 8.4 (19.3) 40 distinct values : 290 381
## [integer] min < med < max: : (43.2%) (56.8%)
## 0 < 3.5 < 192 :
## IQR (CV) : 9 (2.3) :
## : .
##
## 12 magsulf 1. 0 367 (86.6%) IIIIIIIIIIIIIIIII 424 247
## [factor] 2. 1 57 (13.4%) II (63.2%) (36.8%)
##
## 13 meth 1. 0 318 (56.3%) IIIIIIIIIII 565 106
## [factor] 2. 1 247 (43.7%) IIIIIIII (84.2%) (15.8%)
##
## 14 toc 1. 0 438 (77.5%) IIIIIIIIIIIIIII 565 106
## [factor] 2. 1 127 (22.5%) IIII (84.2%) (15.8%)
##
## 15 delivery 1. abdominal 314 (48.4%) IIIIIIIII 649 22
## [factor] 2. vaginal 335 (51.6%) IIIIIIIIII (96.7%) (3.3%)
##
## 16 apg1 Mean (sd) : 4.9 (2.6) 0 : 5 ( 0.8%) 637 34
## [integer] min < med < max: 1 : 92 (14.4%) II (94.9%) (5.1%)
## 0 < 5 < 9 2 : 69 (10.8%) II
## IQR (CV) : 5 (0.5) 3 : 57 ( 8.9%) I
## 4 : 45 ( 7.1%) I
## 5 : 69 (10.8%) II
## 6 : 80 (12.6%) II
## 7 : 80 (12.6%) II
## 8 : 103 (16.2%) III
## 9 : 37 ( 5.8%) I
##
## 17 vent 1. 0 269 (42.0%) IIIIIIII 641 30
## [factor] 2. 1 372 (58.0%) IIIIIIIIIII (95.5%) (4.5%)
##
## 18 pneumo 1. 0 518 (80.3%) IIIIIIIIIIIIIIII 645 26
## [factor] 2. 1 127 (19.7%) III (96.1%) (3.9%)
##
## 19 pda 1. 0 508 (79.1%) IIIIIIIIIIIIIII 642 29
## [factor] 2. 1 134 (20.9%) IIII (95.7%) (4.3%)
##
## 20 cld 1. 0 442 (73.1%) IIIIIIIIIIIIII 605 66
## [factor] 2. 1 163 (26.9%) IIIII (90.2%) (9.8%)
##
## 21 pvh 1. absent 360 (68.4%) IIIIIIIIIIIII 526 145
## [factor] 2. definite 125 (23.8%) IIII (78.4%) (21.6%)
## 3. possible 41 ( 7.8%) I
##
## 22 ivh 1. absent 442 (83.9%) IIIIIIIIIIIIIIII 527 144
## [factor] 2. definite 75 (14.2%) II (78.5%) (21.5%)
## 3. possible 10 ( 1.9%)
##
## 23 ipe 1. absent 472 (89.6%) IIIIIIIIIIIIIIIII 527 144
## [factor] 2. definite 38 ( 7.2%) I (78.5%) (21.5%)
## 3. possible 17 ( 3.2%)
##
## 24 year Mean (sd) : 84.8 (1.6) 542 distinct values : : 650 21
## [numeric] min < med < max: : : : : . : (96.9%) (3.1%)
## 81.5 < 84.9 < 87.5 : : : : : : : : :
## IQR (CV) : 2.6 (0) : : : : : : : : : :
## : : : : : : : : : :
##
## 25 sex 1. female 320 (49.2%) IIIIIIIII 650 21
## [factor] 2. male 330 (50.8%) IIIIIIIIII (96.9%) (3.1%)
##
## 26 dead 1. 0 527 (78.5%) IIIIIIIIIIIIIII 671 0
## [factor] 2. 1 144 (21.5%) IIII (100.0%) (0.0%)
## --------------------------------------------------------------------------------------------------------------Что-то странное происходит в hospstay, предполагаю, что младенцы не умеют путешествовать во времени. В этих (exit, pltct, lol) переменных так же есть экстремальные значения, природа которых мне неизвестна, поэтому удалю их по правилу трех IQR.
NB: после удаления аутлаеров по hospstay исчезли аутлаеры у exit
data <- data %>%
mutate(ID = c(1:nrow(data))) %>%
select(ID, everything()) data <- data %>%
filter(hospstay > 0 | is.na(hospstay),
hospstay < 3 * IQR(hospstay, na.rm = T) | is.na(hospstay),
pltct < 3 * IQR(pltct, na.rm = T) | is.na(pltct),
lol < 3 * IQR(lol, na.rm = T) | is.na(lol))
data %>%
dfSummary()## Data Frame Summary
## Dimensions: 569 x 27
## Duplicates: 0
##
## ---------------------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ----------- ---------------------------- --------------------- ---------------------- ---------- ---------
## 1 ID Mean (sd) : 332.9 (193.8) 569 distinct values . : : . : . . . . . 569 0
## [integer] min < med < max: : : : : : : : : : : (100.0%) (0.0%)
## 1 < 330 < 671 : : : : : : : : : :
## IQR (CV) : 336 (0.6) : : : : : : : : : :
## : : : : : : : : : :
##
## 2 birth Mean (sd) : 84.7 (1.6) 460 distinct values : . 548 21
## [numeric] min < med < max: : : : : . (96.3%) (3.7%)
## 81.5 < 84.9 < 87.5 : : : : : : : : :
## IQR (CV) : 2.5 (0) : : : : : : : : : :
## : : : : : : : : : :
##
## 3 exit Mean (sd) : 84.8 (1.6) 460 distinct values : 538 31
## [numeric] min < med < max: : . : : : (94.6%) (5.4%)
## 81.5 < 84.9 < 87.8 : : : : : : : : .
## IQR (CV) : 2.6 (0) . : : : : : : : : :
## : : : : : : : : : :
##
## 4 hospstay Mean (sd) : 41.6 (30.6) 117 distinct values : 538 31
## [integer] min < med < max: : : : . (94.6%) (5.4%)
## 1 < 36 < 137 : : : : .
## IQR (CV) : 43 (0.7) : : : : : .
## : : : : : : : . . .
##
## 5 lowph Mean (sd) : 7.2 (0.1) 71 distinct values : . 516 53
## [numeric] min < med < max: : : (90.7%) (9.3%)
## 6.5 < 7.2 < 7.5 . : :
## IQR (CV) : 0.2 (0) . : : : :
## . : : : : : .
##
## 6 pltct Mean (sd) : 192.1 (68.1) 222 distinct values : 509 60
## [integer] min < med < max: : : (89.5%) (10.5%)
## 16 < 199 < 324 : : : .
## IQR (CV) : 99 (0.4) . : : : :
## . : : : : : :
##
## 7 race 1. black 303 (55.7%) IIIIIIIIIII 544 25
## [factor] 2. native American 16 ( 2.9%) (95.6%) (4.4%)
## 3. oriental 3 ( 0.6%)
## 4. white 222 (40.8%) IIIIIIII
##
## 8 bwt Mean (sd) : 1110.1 (258.9) 134 distinct values : 567 2
## [integer] min < med < max: . : : : (99.6%) (0.4%)
## 400 < 1150 < 1580 . : : : : :
## IQR (CV) : 400 (0.2) . : : : : : : .
## . : : : : : : : : :
##
## 9 gest Mean (sd) : 29.1 (2.5) 23 distinct values . : 566 3
## [numeric] min < med < max: : : : (99.5%) (0.5%)
## 23 < 29 < 40 : : :
## IQR (CV) : 4 (0.1) : : : :
## . : : : : :
##
## 10 inout 1. born at Duke 462 (81.5%) IIIIIIIIIIIIIIII 567 2
## [factor] 2. transported 105 (18.5%) III (99.6%) (0.4%)
##
## 11 twn 1. 0 432 (78.5%) IIIIIIIIIIIIIII 550 19
## [factor] 2. 1 118 (21.5%) IIII (96.7%) (3.3%)
##
## 12 lol Mean (sd) : 5.2 (6.3) 26 distinct values : 242 327
## [integer] min < med < max: : (42.5%) (57.5%)
## 0 < 3 < 26 :
## IQR (CV) : 8 (1.2) : .
## : : . . .
##
## 13 magsulf 1. 0 306 (85.2%) IIIIIIIIIIIIIIIII 359 210
## [factor] 2. 1 53 (14.8%) II (63.1%) (36.9%)
##
## 14 meth 1. 0 263 (55.4%) IIIIIIIIIII 475 94
## [factor] 2. 1 212 (44.6%) IIIIIIII (83.5%) (16.5%)
##
## 15 toc 1. 0 372 (78.5%) IIIIIIIIIIIIIII 474 95
## [factor] 2. 1 102 (21.5%) IIII (83.3%) (16.7%)
##
## 16 delivery 1. abdominal 274 (50.1%) IIIIIIIIII 547 22
## [factor] 2. vaginal 273 (49.9%) IIIIIIIII (96.1%) (3.9%)
##
## 17 apg1 Mean (sd) : 5 (2.6) 0 : 3 ( 0.6%) 536 33
## [integer] min < med < max: 1 : 74 (13.8%) II (94.2%) (5.8%)
## 0 < 5 < 9 2 : 53 ( 9.9%) I
## IQR (CV) : 4 (0.5) 3 : 46 ( 8.6%) I
## 4 : 35 ( 6.5%) I
## 5 : 61 (11.4%) II
## 6 : 68 (12.7%) II
## 7 : 71 (13.2%) II
## 8 : 96 (17.9%) III
## 9 : 29 ( 5.4%) I
##
## 18 vent 1. 0 236 (43.5%) IIIIIIII 543 26
## [factor] 2. 1 307 (56.5%) IIIIIIIIIII (95.4%) (4.6%)
##
## 19 pneumo 1. 0 441 (81.1%) IIIIIIIIIIIIIIII 544 25
## [factor] 2. 1 103 (18.9%) III (95.6%) (4.4%)
##
## 20 pda 1. 0 441 (81.1%) IIIIIIIIIIIIIIII 544 25
## [factor] 2. 1 103 (18.9%) III (95.6%) (4.4%)
##
## 21 cld 1. 0 393 (76.2%) IIIIIIIIIIIIIII 516 53
## [factor] 2. 1 123 (23.8%) IIII (90.7%) (9.3%)
##
## 22 pvh 1. absent 307 (69.3%) IIIIIIIIIIIII 443 126
## [factor] 2. definite 100 (22.6%) IIII (77.9%) (22.1%)
## 3. possible 36 ( 8.1%) I
##
## 23 ivh 1. absent 378 (85.3%) IIIIIIIIIIIIIIIII 443 126
## [factor] 2. definite 57 (12.9%) II (77.9%) (22.1%)
## 3. possible 8 ( 1.8%)
##
## 24 ipe 1. absent 399 (90.1%) IIIIIIIIIIIIIIIIII 443 126
## [factor] 2. definite 30 ( 6.8%) I (77.9%) (22.1%)
## 3. possible 14 ( 3.2%)
##
## 25 year Mean (sd) : 84.7 (1.6) 460 distinct values : . 548 21
## [numeric] min < med < max: : : : : . (96.3%) (3.7%)
## 81.5 < 84.9 < 87.5 : . : : : : : : :
## IQR (CV) : 2.5 (0) . : : : : : : : : :
## : : : : : : : : : :
##
## 26 sex 1. female 277 (50.5%) IIIIIIIIII 548 21
## [factor] 2. male 271 (49.5%) IIIIIIIII (96.3%) (3.7%)
##
## 27 dead 1. 0 462 (81.2%) IIIIIIIIIIIIIIII 569 0
## [factor] 2. 1 107 (18.8%) III (100.0%) (0.0%)
## ---------------------------------------------------------------------------------------------------------------data %>%
select_if(is.numeric) %>%
select(-ID) %>%
pivot_longer(cols = everything(), names_to = "variable", values_to = "value") %>%
ggplot(aes(x = value)) +
geom_histogram(bins = 10, fill = "skyblue", color = "black") +
facet_wrap(~variable, scales = "free_x") +
theme_minimal() +
theme(axis.text.x = element_text(size=7))## Warning: Removed 582 rows containing non-finite outside the scale range
## (`stat_bin()`).
Теперь выглядит приятнее :)
Task 3
Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?
data %>% group_by(inout) %>% select(inout, lowph) %>% dfSummary()## Data Frame Summary
## Group: inout = born at Duke
## Dimensions: 462 x 2
## Duplicates: 396
##
## -------------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ----------- ----------------------- -------------------- --------------------- --------- ---------
## 2 lowph Mean (sd) : 7.2 (0.1) 65 distinct values : : 420 42
## [numeric] min < med < max: : : (90.9%) (9.1%)
## 6.5 < 7.2 < 7.5 . : : .
## IQR (CV) : 0.2 (0) : : : :
## . : : : : : .
## -------------------------------------------------------------------------------------------------------
##
## Group: inout = transported
## Dimensions: 105 x 2
## Duplicates: 60
##
## ---------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ----------- ----------------------- -------------------- ----------------- --------- ---------
## 2 lowph Mean (sd) : 7.1 (0.2) 44 distinct values : 96 9
## [numeric] min < med < max: : . (91.4%) (8.6%)
## 6.7 < 7.2 < 7.5 : :
## IQR (CV) : 0.2 (0) : : : : .
## . . . : : : : :
## ---------------------------------------------------------------------------------------------------
##
## Group: inout = NA
## Dimensions: 2 x 2
## Duplicates: 1
##
## ----------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ----------- ---------------- -------------------- ------- -------- ----------
## 2 lowph All NA's 0 2
## [numeric] (0.0%) (100.0%)
## ----------------------------------------------------------------------------------data %>%
t_test(lowph ~ inout, detailed = T)## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "born at Duke" - "transported", or divided in the order "born at Duke" /
## "transported" for ratio-based statistics. To specify this order yourself,
## supply `order = c("born at Duke", "transported")`.## # A tibble: 1 × 7
## statistic t_df p_value alternative estimate lower_ci upper_ci
## <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 5.20 130. 0.000000744 two.sided 0.0881 0.0546 0.122У транспортированных пациентов значение рН меньше, чем у рожденных в Дюке (MD=0.088, CI 95%: [0.055; 0.122], p < \(10^{-6}\)). Поскольку мы знаем, что более низкое значение рН ассоциировано с более высокой смертностью, можно предположить, что транспортировка пациентов может быть также связана с более высокой смертностью.
Task 4
Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.
corr_data <- data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
psych::corr.test(method = "spearman", adjust="BH")
cor_melted <- reshape2::melt(corr_data$r)
p_melted <- reshape2::melt(corr_data$p)
cor_plot_data <- merge(cor_melted, p_melted, by = c("Var1", "Var2"))
colnames(cor_plot_data) <- c("Var1", "Var2", "Correlation", "Adjusted_P")
ggplot(cor_plot_data, aes(Var1, Var2)) +
geom_tile(aes(fill = Correlation), color = "white") +
scale_fill_gradient2(low = "#433CE1", high = "#E13C59", mid = "grey90",
midpoint = 0, limit = c(-1, 1), name = "Correlation") +
geom_text(aes(label = ifelse(Adjusted_P < 0.05, sprintf("%.2f", Adjusted_P), "")),
color = "black", size = 4) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
axis.text.y = element_text(size = 12)) +
labs(x = "", y = "")
data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
ggpairs(progress = FALSE) +
theme_bw() +
theme(axis.text = element_text(size=6))
Task 5
Постройте иерархическую кластеризацию на этом датафрейме.
rownames(data) <- data$ID
res <- data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
na.omit() %>%
scale() %>%
dist(method = "euclidean") %>%
hclust(method = "ward.D2")
res %>%
fviz_dend(cex = 0.2,
k = 3,
k_colors = "jco")## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Task 6
Сделайте одновременный график heatmap и иерархической кластеризации. Интерпретируйте результат.
data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
na.omit() %>%
scale() %>%
pheatmap::pheatmap(color=colorRampPalette(c("#433CE1", "grey90", "#E13C59"))(50),
fontsize_row = 0.00001, fontsize_col = 7, angle_col = 0)
data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
na.omit() %>%
scale() %>%
pheatmap::pheatmap(color=colorRampPalette(c("#433CE1", "grey90", "#E13C59"))(50),
fontsize_row = 0.00001, fontsize_col = 7, angle_col = 0, kmeans_k=4)
Визульно, можно выделить на хитмапе 1 4 кластера, чтобы их проинтерпретировать, можно усреднить значения в кластере, как на хитмапе 2. Первые 2 кластера довольно сильно похожи друг на друга, за исключением того, что пациенты из 1 кластера имели более выскойи уровень тромбоцитов и скор APG в отличие от пациентов из кластера 1. Пациенты из кластера 3 в отличие от остальных имеют наиболее длительный срок нахождения в больнице. Отличительными чертами пациентов из кластера 4 является самый низкий уровень рН, масса тела при рождении, а также минимальные сроки беременности.
Task 7 & 8
Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?
Постройте biplot график для PCA. Раскрасьте его по значению колонки ‘dead’.
pca.data <- data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
na.omit() %>%
scale() %>%
prcomp()
fviz_eig(pca.data,
addlabels = T,
ylim = c(0, 50))
В нашем случае необходимо проводить шкалирование данных, так как все переменные имеют разные единицы измерения и разный размах. Первые две компоненты довольно хорошо объясняют данные, объясняя около 55% вариабельности.
pca.data$x %>%
ggplot() +
geom_point(aes(x = PC1, y = PC2)) +
theme_bw()
Вклад каждой переменной в компоненту:
pcs <- reshape2::melt(pca.data$rotation)
ggplot(pcs, aes(Var1, Var2)) +
geom_tile(aes(fill = value), color = "white") +
scale_fill_gradient2(low = "#433CE1", high = "#E13C59", mid = "grey90",
midpoint = 0, limit = c(-1, 1)) +
geom_text(aes(label = sprintf("%.2f", value)),
color = "black", size = 4) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
axis.text.y = element_text(size = 12)) +
labs(x = "", y = "")
В первую компоненту основной положительный вклад имеет время,
проведенное в госпитале, и негативная масса тела при рождении и время
беременности (что похоже на последний кластер в анализе выше).
Во вторую компоненту основной вклад вносит продолжительность родов. В третий уровень тромбоцитов и APG скор в большей степени и в меньшей степени продолжительность нахождения в больнице.
deaths <- data %>% filter(ID %in% rownames(pca.data$x)) %>% select(dead)
levels(deaths$dead) <- c("Alive", "Dead")
ggbiplot(pca.data, choices=1:2,
scale=1, alpha = 0.5, groups = as.factor(deaths$dead)) +
theme_bw()
Task 9
Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.
explained_variance <- pca.data$sdev
explained_variance_pct <- round(100 * explained_variance ** 2 / sum(explained_variance ** 2), 2)
comp <- pca.data$rotation
loadings <- comp
for (i in seq(explained_variance)){
loadings[,i] <- comp[,i] * explained_variance[i]
}
features <- rownames(pca.data$rotation)
fig <- plot_ly(
data = as.data.frame(pca.data$x), x = ~PC1, y = ~PC2,
text = ~paste("ID: ", rownames(pca.data$x)),
color = ~as.factor(deaths$dead),
colors = c('#636EFA','#EF553B'),
type = 'scatter',
mode = 'markers') %>%
layout(
legend=list(title=list(text='color')),
plot_bgcolor = "#e5ecf6",
xaxis = list(
title = sprintf("standardized PC1 (%.2f%% explained var.)", explained_variance_pct[1])),
yaxis = list(
title = sprintf("standardized PC2 (%.2f%% explained var.)", explained_variance_pct[2])))
for (i in seq(length(features))){
fig <- fig %>%
add_segments(x = 0,
xend = loadings[i, 1],
y = 0,
yend = loadings[i, 2],
line = list(color = 'black'),
inherit = FALSE,
showlegend = FALSE) %>%
add_annotations(x=loadings[i, 1],
y=loadings[i, 2],
ax = 0,
ay = 0,
text = features[i],
xanchor = 'center',
yanchor= 'bottom')
}
figTask 10
Дайте содержательную интерпретацию PCA анализу. Почему использовать колонку ‘dead’ для выводов об ассоциации с выживаемостью некорректно?
To be done
Task 11
Приведите ваши данные к размерности в две колонки через UMAP. Сравните результаты отображения точек между алгоритмами PCA и UMAP
umap_comps <- data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
na.omit() %>%
recipe(~.) %>%
step_normalize(all_predictors()) %>%
step_umap(all_predictors()) %>%
prep() %>%
juice()
combined_data <- cbind(umap_comps, pca.data$x[,1:2], deaths)umap_plot <- combined_data %>%
ggplot(aes(UMAP1, UMAP2)) +
geom_point(aes(color = as.character(dead)),
alpha = 0.7, size = 2) +
labs(color = NULL) +
theme_bw()
pca_plot <- combined_data %>%
ggplot(aes(PC1, PC2)) +
geom_point(aes(color = as.character(dead)),
alpha = 0.7, size = 2) +
labs(color = NULL) +
theme_bw()
ggarrange(
pca_plot, umap_plot,
nrow = 1, ncol = 2,
widths = c(1, 1)
)
В целом, результаты получились похожими, хотя умершие субъекты лучше
“кластеризуются” в UMAP.